sn-ENS neuron integration

advanced:

clustering modification

cell composition

all markers

check signature score of some genesets

source("/Shared_win/projects/RNA_normal/analysis.10x.r")

load integrated obj

GEX.seur <- readRDS("./sn10x_WYS.sct_anno.rds")
GEX.seur
## An object of class Seurat 
## 47356 features across 19649 samples within 3 assays 
## Active assay: integrated (2397 features, 2397 variable features)
##  2 other assays present: RNA, SCT
##  3 dimensional reductions calculated: pca, tsne, umap
#scales::show_col(ggsci::pal_igv("default")(49))
color.cnt <- scales::hue_pal()(4)[c(2,1,4,3)]
color.cnt
## [1] "#7CAE00" "#F8766D" "#C77CFF" "#00BFC4"
head(GEX.seur@meta.data)
##                        orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1_1 Nb5d.PBS_INF       3257         1801 0.36843721  0.3991403
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS_INF       1511          997 0.66181337  0.4632694
## AAACCCAGTTTGTTCT-1_1 Nb5d.PBS_INF       2855         1577 0.98073555  0.3152364
## AAACCCATCCTAGCCT-1_1 Nb5d.PBS_INF       2433         1451 0.08220304  0.3699137
## AAACCCATCGAAACAA-1_1 Nb5d.PBS_INF       3129         1656 0.12783637  0.4474273
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS_INF       2201         1294 0.22716947  0.2271695
##                           S.Score     G2M.Score Phase      cnt  rep newAnno
## AAACCCACAAGACGAC-1_1  0.011590405 -0.0004169865     S Nb5d.INF rep4    EMN3
## AAACCCAGTGGGCTCT-1_1 -0.024203070  0.0012414826   G2M Nb5d.PBS rep4   IPAN1
## AAACCCAGTTTGTTCT-1_1 -0.013980476  0.0039329410   G2M Nb5d.INF rep1    EMN3
## AAACCCATCCTAGCCT-1_1 -0.028925620 -0.0132582758    G1 Nb5d.INF rep2    EMN1
## AAACCCATCGAAACAA-1_1 -0.008077289 -0.0028336129    G1 Nb5d.PBS rep3   IPAN4
## AAACCCATCGGTCAGC-1_1 -0.023612751  0.0327239644   G2M Nb5d.PBS rep4    EMN1
##                         sample tissue nCount_SCT nFeature_SCT condition
## AAACCCACAAGACGAC-1_1 Nb5d.INF4  Ileum       2598         1771   INF_CTL
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS4  Ileum       1696          971   PBS_CTL
## AAACCCAGTTTGTTCT-1_1 Nb5d.INF1  Ileum       2493         1558   INF_CTL
## AAACCCATCCTAGCCT-1_1 Nb5d.INF2  Ileum       2319         1427   INF_CTL
## AAACCCATCGAAACAA-1_1 Nb5d.PBS3  Ileum       2543         1620   PBS_CTL
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS4  Ileum       2166         1274   PBS_CTL
##                      integrated_snn_res.2 seurat_clusters
## AAACCCACAAGACGAC-1_1                    8               8
## AAACCCAGTGGGCTCT-1_1                    7               7
## AAACCCAGTTTGTTCT-1_1                    8               8
## AAACCCATCCTAGCCT-1_1                    2               2
## AAACCCATCGAAACAA-1_1                   18              18
## AAACCCATCGGTCAGC-1_1                    2               2
DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 3.25, group.by = "newAnno") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 3.75, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)

DimPlot(subset(GEX.seur, subset = cnt %in% c("Nb5d.PBS","Nb5d.INF")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[1:2]) +
DimPlot(subset(GEX.seur, subset = cnt %in% c("Nb5d.CTL","Nb5d.CKO")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[3:4])

re-clustering

DefaultAssay(GEX.seur) <- "integrated"
GEX.seur
## An object of class Seurat 
## 47356 features across 19649 samples within 3 assays 
## Active assay: integrated (2397 features, 2397 variable features)
##  2 other assays present: RNA, SCT
##  3 dimensional reductions calculated: pca, tsne, umap
ElbowPlot(GEX.seur, ndims = 100)

ElbowPlot(GEX.seur, ndims = 50)

PCsct <- 1:20

GEX.seur <- FindNeighbors(GEX.seur, k.param = 18, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution =2, reduction = 'pca', verbose = T)
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 19649
## Number of edges: 584347
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8315
## Number of communities: 31
## Elapsed time: 2 seconds
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")

DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "newAnno") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")

DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")

DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)

DefaultAssay(GEX.seur) <- "SCT"
GEX.seur
## An object of class Seurat 
## 47356 features across 19649 samples within 3 assays 
## Active assay: SCT (20434 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
# sort_clusters                         

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(0,1,11,5,16,19,        # EMN1
                                            15,8,                  # EMN2
                                            4,10,                  # EMN3
                                            17,                    # EMN4
                                            25,                    # EMN5
                                            6,3,2,20,              # IMN1
                                            14,                    # IMN2
                                            9,28,                  # IMN3
                                            26,                    # IMN4
                                            12,                    # IN1
                                            24,                    # IN2
                                            30,                    # IN3
                                            21,22,13, 29,          # IPAN1
                                            7, 18,                 # IPAN2
                                            27,                    # IPAN3
                                            23                     # IPAN4
                                            ))








# ttAnno1 (based on individual newAnno)             
#    hold 
# ttAnno2 (sub-cluster IPAN1/2 as .1 & .2)                                
#    hold
markers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Itga6","Cntnap5b",
                       "Pgm5","Chgb",
                       "Nxph1",
                       "Gabrb1","Glp2r","Nebl",
                       "Lrrc7",
                       "Ryr3","Eda",
                       "Hgf","Lama2","Efnb2",
                       "Tac1",
                       "Kctd8",
                       "Ptn",
                       "Ntrk2","Penk","Itgb8",
                       "Fut9","Nfatc1","Egfr","Pdia5",
                       "Ahr","Mgll",
                       "Aff3",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2",
                       "Col25a1",
                       "Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5","Mid1","Igfbp5",
                       "Ppara",
                       "Pcdh11x","Adcy8","Grp"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Nova1",
                      "Cdh10","Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Vipr2","Adamts9","Sst","Kcnn2"
                      ),
                 IPAN=c("Calb2","Adcy1","Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Itgb6","Met","Itgbl1",
                        
                        "Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9",
                        "Car10","Dcc",
                        "Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))
check.markers.ss <- as.vector(unlist(markers.new.ss))
check.markers.ss
##   [1] "Chat"     "Bnc2"     "Tox"      "Ptprt"    "Gfra2"    "Oprk1"   
##   [7] "Adamtsl1" "Fbxw15"   "Fbxw24"   "Chrna7"   "Satb1"    "Itga6"   
##  [13] "Cntnap5b" "Pgm5"     "Chgb"     "Nxph1"    "Gabrb1"   "Glp2r"   
##  [19] "Nebl"     "Lrrc7"    "Ryr3"     "Eda"      "Hgf"      "Lama2"   
##  [25] "Efnb2"    "Tac1"     "Kctd8"    "Ptn"      "Ntrk2"    "Penk"    
##  [31] "Itgb8"    "Fut9"     "Nfatc1"   "Egfr"     "Pdia5"    "Ahr"     
##  [37] "Mgll"     "Aff3"     "Chrm3"    "Nos1"     "Kcnab1"   "Gfra1"   
##  [43] "Etv1"     "Man1a"    "Airn"     "Adcy2"    "Col25a1"  "Cmah"    
##  [49] "Creb5"    "Vip"      "Pde1a"    "Ebf1"     "Gpc5"     "Mid1"    
##  [55] "Igfbp5"   "Ppara"    "Pcdh11x"  "Adcy8"    "Grp"      "Npas3"   
##  [61] "Synpr"    "St18"     "Gal"      "Nova1"    "Cdh10"    "Kcnk13"  
##  [67] "Neurod6"  "Moxd1"    "Sctr"     "Piezo1"   "Vipr2"    "Adamts9" 
##  [73] "Sst"      "Kcnn2"    "Calb2"    "Adcy1"    "Calcb"    "Nmu"     
##  [79] "Adgrg6"   "Pcdh10"   "Ngfr"     "Galr1"    "Il7"      "Aff2"    
##  [85] "Gpr149"   "Itgb6"    "Met"      "Itgbl1"   "Cdh6"     "Cdh8"    
##  [91] "Clstn2"   "Ano2"     "Ntrk3"    "Cpne4"    "Vwc2l"    "Cdh9"    
##  [97] "Car10"    "Dcc"      "Scgn"     "Vcan"     "Cck"      "Piezo2"  
## [103] "Kcnh7"    "Rerg"     "Bmpr1b"   "Skap1"    "Ntng1"    "Tafa2"   
## [109] "Nxph2"    "Actl6b"   "Phox2b"   "Grid2"    "Ide"      "Btaf1"   
## [115] "Slc15a2"  "Ccser1"   "Hdac9"    "Rspo2"    "Grem2"
length(check.markers.ss)
## [1] 119
sum(check.markers.ss %in% rownames(GEX.seur))
## [1] 119
pm.All_new.a <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "sort_clusters",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - sort_clusters")
pm.All_new.a

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "sort_clusters")

mod1

DefaultAssay(GEX.seur) <- "integrated"
GEX.seur
## An object of class Seurat 
## 47356 features across 19649 samples within 3 assays 
## Active assay: integrated (2397 features, 2397 variable features)
##  2 other assays present: RNA, SCT
##  3 dimensional reductions calculated: pca, tsne, umap
# remove C29
GEX.seur <- subset(GEX.seur, subset= seurat_clusters %in% c(0:28,30))
GEX.seur
## An object of class Seurat 
## 47356 features across 19460 samples within 3 assays 
## Active assay: integrated (2397 features, 2397 variable features)
##  2 other assays present: RNA, SCT
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")

ElbowPlot(GEX.seur, ndims = 50)

PCsct <- 1:18

GEX.seur <- FindNeighbors(GEX.seur, k.param = 12, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution =2.5, reduction = 'pca', verbose = T)
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 19460
## Number of edges: 442022
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8093
## Number of communities: 35
## Elapsed time: 2 seconds
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "seurat_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "seurat_clusters")

mod2

# remove C34
GEX.seur <- subset(GEX.seur, subset= seurat_clusters %in% c(0:33))
GEX.seur
## An object of class Seurat 
## 47356 features across 19418 samples within 3 assays 
## Active assay: integrated (2397 features, 2397 variable features)
##  2 other assays present: RNA, SCT
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")

ElbowPlot(GEX.seur, ndims = 50)

PCsct <- 1:18

GEX.seur <- FindNeighbors(GEX.seur, k.param = 12, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution =2.5, reduction = 'pca', verbose = T)
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 19418
## Number of edges: 441598
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8095
## Number of communities: 34
## Elapsed time: 2 seconds
pp.umap1 <- DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")
pp.umap1

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "seurat_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "seurat_clusters")

DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "newAnno") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")

DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")

DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)

intAnno

set modified annotations for integration as intAnno1/2

DefaultAssay(GEX.seur) <- "SCT"
GEX.seur
## An object of class Seurat 
## 47356 features across 19418 samples within 3 assays 
## Active assay: SCT (20434 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.seur@meta.data[,grep("snn",colnames(GEX.seur@meta.data))] <- NULL
head(GEX.seur@meta.data)
##                        orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1_1 Nb5d.PBS_INF       3257         1801 0.36843721  0.3991403
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS_INF       1511          997 0.66181337  0.4632694
## AAACCCAGTTTGTTCT-1_1 Nb5d.PBS_INF       2855         1577 0.98073555  0.3152364
## AAACCCATCCTAGCCT-1_1 Nb5d.PBS_INF       2433         1451 0.08220304  0.3699137
## AAACCCATCGAAACAA-1_1 Nb5d.PBS_INF       3129         1656 0.12783637  0.4474273
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS_INF       2201         1294 0.22716947  0.2271695
##                           S.Score     G2M.Score Phase      cnt  rep newAnno
## AAACCCACAAGACGAC-1_1  0.011590405 -0.0004169865     S Nb5d.INF rep4    EMN3
## AAACCCAGTGGGCTCT-1_1 -0.024203070  0.0012414826   G2M Nb5d.PBS rep4   IPAN1
## AAACCCAGTTTGTTCT-1_1 -0.013980476  0.0039329410   G2M Nb5d.INF rep1    EMN3
## AAACCCATCCTAGCCT-1_1 -0.028925620 -0.0132582758    G1 Nb5d.INF rep2    EMN1
## AAACCCATCGAAACAA-1_1 -0.008077289 -0.0028336129    G1 Nb5d.PBS rep3   IPAN4
## AAACCCATCGGTCAGC-1_1 -0.023612751  0.0327239644   G2M Nb5d.PBS rep4    EMN1
##                         sample tissue nCount_SCT nFeature_SCT condition
## AAACCCACAAGACGAC-1_1 Nb5d.INF4  Ileum       2598         1771   INF_CTL
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS4  Ileum       1696          971   PBS_CTL
## AAACCCAGTTTGTTCT-1_1 Nb5d.INF1  Ileum       2493         1558   INF_CTL
## AAACCCATCCTAGCCT-1_1 Nb5d.INF2  Ileum       2319         1427   INF_CTL
## AAACCCATCGAAACAA-1_1 Nb5d.PBS3  Ileum       2543         1620   PBS_CTL
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS4  Ileum       2166         1274   PBS_CTL
##                      seurat_clusters sort_clusters
## AAACCCACAAGACGAC-1_1              11            10
## AAACCCAGTGGGCTCT-1_1              22            21
## AAACCCAGTTTGTTCT-1_1              11             8
## AAACCCATCCTAGCCT-1_1               4             1
## AAACCCATCGAAACAA-1_1              19            23
## AAACCCATCGGTCAGC-1_1               8             1
# sort_clusters                         

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(0,5,4,8,16,17,14,        # EMN1
                                            10,27,11,               # EMN2
                                            2,6,                    # EMN3
                                            28,                     # EMN4
                                            26,                     # EMN5
                                            1,9,15,3,18,              # IMN1
                                            13,24,                    # IMN2
                                            20,                    # IMN3
                                            29,                    # IMN4
                                            12,                    # IN1
                                            25,                    # IN2
                                            33,                    # IN3
                                            22,31, 23,30,          # IPAN1
                                            7, 21,                 # IPAN2
                                            32,                    # IPAN3
                                            19                     # IPAN4
                                            ))








# intAnno1 (based on individual newAnno)             

GEX.seur$intAnno1 <- as.character(GEX.seur$seurat_clusters)

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(0,5,4,8,16,17,14)] <- "EMN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(10,27,11)] <- "EMN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(2,6)] <- "EMN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(28)] <- "EMN4"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(26)] <- "EMN5"

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(1,9,15,3,18)] <- "IMN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(13,24)] <- "IMN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(20)] <- "IMN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(29)] <- "IMN4"

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(12)] <- "IN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(25)] <- "IN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(33)] <- "IN3"

GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(22,31,23,30)] <- "IPAN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(7,21)] <- "IPAN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(32)] <- "IPAN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(19)] <- "IPAN4"

GEX.seur$intAnno1 <- factor(GEX.seur$intAnno1,
                            levels = c(paste0("EMN",1:5),
                                       paste0("IMN",1:4),
                                       paste0("IN",1:3),
                                       paste0("IPAN",1:4)))


# intAnno2 (sub-cluster IPAN1/2 as .1 & .2)                                

GEX.seur$intAnno2 <- as.character(GEX.seur$seurat_clusters)

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(0,5,4,8,16,17,14)] <- "EMN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(10,27,11)] <- "EMN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(2,6)] <- "EMN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(28)] <- "EMN4"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(26)] <- "EMN5"

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(1,9,15,3,18)] <- "IMN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(13,24)] <- "IMN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(20)] <- "IMN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(29)] <- "IMN4"

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(12)] <- "IN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(25)] <- "IN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(33)] <- "IN3"

GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(22,31)] <- "IPAN1.1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(23,30)] <- "IPAN1.2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(7)] <- "IPAN2.1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(21)] <- "IPAN2.2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(32)] <- "IPAN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(19)] <- "IPAN4"

GEX.seur$intAnno2 <- factor(GEX.seur$intAnno2,
                            levels = c(paste0("EMN",1:5),
                                       paste0("IMN",1:4),
                                       paste0("IN",1:3),
                                       paste0("IPAN",c(1.1,1.2,2.1,2.2,3,4))))
pp.umap2 <- DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "intAnno1") +
DimPlot(GEX.seur, label = T, label.size = 2.65, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "intAnno2")
pp.umap2 

pp.umap3 <- DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)
pp.umap3

markers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
                       "Gfra2","Oprk1","Adamtsl1", 
                       "Fbxw15","Fbxw24","Chrna7",
                       "Satb1","Itga6","Cntnap5b",
                       "Pgm5","Chgb",
                       "Nxph1",
                       "Gabrb1","Glp2r","Nebl",
                       "Lrrc7",
                       "Ryr3","Eda",
                       "Hgf","Lama2","Efnb2",
                       "Tac1",
                       "Kctd8",
                       "Ptn",
                       "Ntrk2","Penk","Itgb8",
                       "Fut9","Nfatc1","Egfr","Pdia5",
                       "Ahr","Mgll",
                       "Aff3",
                       "Chrm3"
                       ),
                 IMN=c("Nos1","Kcnab1",
                       "Gfra1","Etv1",
                       "Man1a","Airn",
                       "Adcy2",
                       "Col25a1",
                       "Cmah","Creb5","Vip","Pde1a",
                       "Ebf1","Gpc5","Mid1","Igfbp5",
                       "Ppara",
                       "Pcdh11x","Adcy8","Grp"
                       ),
                 IN=c("Npas3","Synpr","St18","Gal",
                      "Nova1",
                      "Cdh10","Kcnk13",
                      "Neurod6","Moxd1","Sctr",
                      "Piezo1","Vipr2","Adamts9","Sst","Kcnn2"
                      ),
                 IPAN=c("Calb2","Adcy1","Calcb","Nmu","Adgrg6","Pcdh10",
                        "Ngfr","Galr1","Il7","Aff2",
                        "Gpr149","Itgb6","Met","Itgbl1",
                        
                        "Cdh6","Cdh8",
                        "Clstn2","Ano2","Ntrk3",
                        "Cpne4","Vwc2l","Cdh9",
                        "Car10","Dcc",
                        "Scgn",
                        "Vcan","Cck","Piezo2","Kcnh7",
                        "Rerg","Bmpr1b","Skap1","Ntng1",
                        "Tafa2","Nxph2"),
                 CONTAM=c("Actl6b","Phox2b"),
                 pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))
check.markers.ss <- as.vector(unlist(markers.new.ss))
check.markers.ss
##   [1] "Chat"     "Bnc2"     "Tox"      "Ptprt"    "Gfra2"    "Oprk1"   
##   [7] "Adamtsl1" "Fbxw15"   "Fbxw24"   "Chrna7"   "Satb1"    "Itga6"   
##  [13] "Cntnap5b" "Pgm5"     "Chgb"     "Nxph1"    "Gabrb1"   "Glp2r"   
##  [19] "Nebl"     "Lrrc7"    "Ryr3"     "Eda"      "Hgf"      "Lama2"   
##  [25] "Efnb2"    "Tac1"     "Kctd8"    "Ptn"      "Ntrk2"    "Penk"    
##  [31] "Itgb8"    "Fut9"     "Nfatc1"   "Egfr"     "Pdia5"    "Ahr"     
##  [37] "Mgll"     "Aff3"     "Chrm3"    "Nos1"     "Kcnab1"   "Gfra1"   
##  [43] "Etv1"     "Man1a"    "Airn"     "Adcy2"    "Col25a1"  "Cmah"    
##  [49] "Creb5"    "Vip"      "Pde1a"    "Ebf1"     "Gpc5"     "Mid1"    
##  [55] "Igfbp5"   "Ppara"    "Pcdh11x"  "Adcy8"    "Grp"      "Npas3"   
##  [61] "Synpr"    "St18"     "Gal"      "Nova1"    "Cdh10"    "Kcnk13"  
##  [67] "Neurod6"  "Moxd1"    "Sctr"     "Piezo1"   "Vipr2"    "Adamts9" 
##  [73] "Sst"      "Kcnn2"    "Calb2"    "Adcy1"    "Calcb"    "Nmu"     
##  [79] "Adgrg6"   "Pcdh10"   "Ngfr"     "Galr1"    "Il7"      "Aff2"    
##  [85] "Gpr149"   "Itgb6"    "Met"      "Itgbl1"   "Cdh6"     "Cdh8"    
##  [91] "Clstn2"   "Ano2"     "Ntrk3"    "Cpne4"    "Vwc2l"    "Cdh9"    
##  [97] "Car10"    "Dcc"      "Scgn"     "Vcan"     "Cck"      "Piezo2"  
## [103] "Kcnh7"    "Rerg"     "Bmpr1b"   "Skap1"    "Ntng1"    "Tafa2"   
## [109] "Nxph2"    "Actl6b"   "Phox2b"   "Grid2"    "Ide"      "Btaf1"   
## [115] "Slc15a2"  "Ccser1"   "Hdac9"    "Rspo2"    "Grem2"
length(check.markers.ss)
## [1] 119
sum(check.markers.ss %in% rownames(GEX.seur))
## [1] 119
pm.All_new.a <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "sort_clusters",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - sort_clusters")
pm.All_new.a

pm.All_new.b <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "intAnno1",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno1")
pm.All_new.b

pm.All_new.c <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno2")
pm.All_new.c

#
pm.All_new.c1 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.PBS")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.PBS - intAnno2")
pm.All_new.c1

#
pm.All_new.c2 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.INF")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.INF - intAnno2")
pm.All_new.c2

#
pm.All_new.c3 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.CTL")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.CTL - intAnno2")
pm.All_new.c3

#
pm.All_new.c4 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.CKO")), features = check.markers.ss, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.CKO - intAnno2")
pm.All_new.c4

#saveRDS(GEX.seur, "./sn10x_WYS.sct_anno.s.rds")
#GEX.seur <- readRDS("./sn10x_WYS.sct_anno.s.rds")

cell composition

pumap.test0 <- DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno2") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pumap.test0

pumap.test0s <- DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pumap.test0s

PBS_INF

test1.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.PBS","Nb5d.INF"))
test1.seur
## An object of class Seurat 
## 47356 features across 8232 samples within 3 assays 
## Active assay: SCT (20434 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
color.test1 <- color.cnt[1:2]
pumap.test1 <- DimPlot(test1.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno2") +
  DimPlot(test1.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test1)
pumap.test1

pumap.test1s <- DimPlot(test1.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno1") +
  DimPlot(test1.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test1)
pumap.test1s

heatmap

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test1.seur$sample,
      clusters=test1.seur$intAnno1)[c(5:8,1:4),],
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test1.seur$sample,
                                clusters=test1.seur$intAnno1)[c(5:8,1:4),] ),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test1.seur$sample,
      clusters=test1.seur$intAnno2)[c(5:8,1:4),],
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test1.seur$sample,
                                clusters=test1.seur$intAnno2)[c(5:8,1:4),] ),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot

stat1_intAnno1 <- test1.seur@meta.data[,c("cnt","rep","intAnno1")]

stat1_intAnno1.s <- stat1_intAnno1 %>%
  group_by(cnt,rep,intAnno1) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom1.intAnno1 <- stat1_intAnno1.s %>%
  ggplot(aes(x = intAnno1, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title="Cell Composition of Nb5d.PBS_INF, intAnno1", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno1, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom1.intAnno1

stat1_intAnno2 <- test1.seur@meta.data[,c("cnt","rep","intAnno2")]

stat1_intAnno2.s <- stat1_intAnno2 %>%
  group_by(cnt,rep,intAnno2) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom1.intAnno2 <- stat1_intAnno2.s %>%
  ggplot(aes(x = intAnno2, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title="Cell Composition of Nb5d.PBS_INF, intAnno2", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno2, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom1.intAnno2

sig.test

glm.nb - anova.LRT

Sort.N <- c("Nb5d.PBS","Nb5d.INF")

glm.nb_anova.1.intAnno1 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat1_intAnno1.s_N <- stat1_intAnno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat1_intAnno1.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat1_intAnno1.s_N$total <- total.N[paste0(stat1_intAnno1.s_N$cnt,
                                            "_",
                                            stat1_intAnno1.s_N$rep),"total"]
      
      glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat1_intAnno1.s_N$intAnno1)){
        glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.1.intAnno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.intAnno1)))
rownames(glm.nb_anova.1.intAnno1_df) <- names(glm.nb_anova.1.intAnno1)
colnames(glm.nb_anova.1.intAnno1_df) <- gsub("X","C",colnames(glm.nb_anova.1.intAnno1_df))
glm.nb_anova.1.intAnno1_df
##                          EMN1      EMN2      EMN3      EMN4      EMN5
## Nb5d.PBSvsNb5d.INF 0.02915675 0.8569496 0.2429792 0.7228933 0.9970758
##                          IMN1      IMN2      IMN3       IMN4      IN1       IN2
## Nb5d.PBSvsNb5d.INF 0.01061606 0.8302173 0.9288441 0.09847105 0.669224 0.2176931
##                          IN3     IPAN1     IPAN2     IPAN3     IPAN4
## Nb5d.PBSvsNb5d.INF 0.6229168 0.4983735 0.4164602 0.7004949 0.6806156
#options (scipen = 999)
round(glm.nb_anova.1.intAnno1_df,4)
##                      EMN1   EMN2  EMN3   EMN4   EMN5   IMN1   IMN2   IMN3
## Nb5d.PBSvsNb5d.INF 0.0292 0.8569 0.243 0.7229 0.9971 0.0106 0.8302 0.9288
##                      IMN4    IN1    IN2    IN3  IPAN1  IPAN2  IPAN3  IPAN4
## Nb5d.PBSvsNb5d.INF 0.0985 0.6692 0.2177 0.6229 0.4984 0.4165 0.7005 0.6806
Sort.N <- c("Nb5d.PBS","Nb5d.INF")

glm.nb_anova.1.intAnno2 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat1_intAnno2.s_N <- stat1_intAnno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat1_intAnno2.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat1_intAnno2.s_N$total <- total.N[paste0(stat1_intAnno2.s_N$cnt,
                                            "_",
                                            stat1_intAnno2.s_N$rep),"total"]
      
      glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat1_intAnno2.s_N$intAnno2)){
        glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.1.intAnno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.intAnno2)))
rownames(glm.nb_anova.1.intAnno2_df) <- names(glm.nb_anova.1.intAnno2)
colnames(glm.nb_anova.1.intAnno2_df) <- gsub("X","C",colnames(glm.nb_anova.1.intAnno2_df))
glm.nb_anova.1.intAnno2_df
##                          EMN1      EMN2      EMN3      EMN4      EMN5
## Nb5d.PBSvsNb5d.INF 0.02915675 0.8569496 0.2429792 0.7228933 0.9970758
##                          IMN1      IMN2      IMN3       IMN4      IN1       IN2
## Nb5d.PBSvsNb5d.INF 0.01061606 0.8302173 0.9288441 0.09847105 0.669224 0.2176931
##                          IN3   IPAN1.1  IPAN1.2   IPAN2.1   IPAN2.2     IPAN3
## Nb5d.PBSvsNb5d.INF 0.6229168 0.8709588 0.400638 0.7016721 0.3976744 0.7004949
##                        IPAN4
## Nb5d.PBSvsNb5d.INF 0.6806156
#options (scipen = 999)
round(glm.nb_anova.1.intAnno2_df,4)
##                      EMN1   EMN2  EMN3   EMN4   EMN5   IMN1   IMN2   IMN3
## Nb5d.PBSvsNb5d.INF 0.0292 0.8569 0.243 0.7229 0.9971 0.0106 0.8302 0.9288
##                      IMN4    IN1    IN2    IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2
## Nb5d.PBSvsNb5d.INF 0.0985 0.6692 0.2177 0.6229   0.871  0.4006  0.7017  0.3977
##                     IPAN3  IPAN4
## Nb5d.PBSvsNb5d.INF 0.7005 0.6806

CTL_CKO

test2.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.CTL","Nb5d.CKO"))
test2.seur
## An object of class Seurat 
## 47356 features across 11186 samples within 3 assays 
## Active assay: SCT (20434 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap
color.test2 <- color.cnt[3:4]
pumap.test2 <- DimPlot(test2.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno2") +
  DimPlot(test2.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test2)
pumap.test2

pumap.test2s <- DimPlot(test2.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno1") +
  DimPlot(test2.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test2)
pumap.test2s

heatmap

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test2.seur$sample,
      clusters=test2.seur$intAnno1)[c(4:7,1:3),],
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test2.seur$sample,
                                clusters=test2.seur$intAnno1)[c(4:7,1:3),] ),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test2.seur$sample,
      clusters=test2.seur$intAnno2)[c(4:7,1:3),],
                   main = "Cell Count",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(cnt=test2.seur$sample,
                                clusters=test2.seur$intAnno2)[c(4:7,1:3),] ),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      gaps_col = c(5,9,12),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

barplot

stat2_intAnno1 <- test2.seur@meta.data[,c("cnt","rep","intAnno1")]

stat2_intAnno1.s <- stat2_intAnno1 %>%
  group_by(cnt,rep,intAnno1) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom2.intAnno1 <- stat2_intAnno1.s %>%
  ggplot(aes(x = intAnno1, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title="Cell Composition of Nb5d.CTL_CKO, intAnno1", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno1, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom2.intAnno1

stat2_intAnno2 <- test2.seur@meta.data[,c("cnt","rep","intAnno2")]

stat2_intAnno2.s <- stat2_intAnno2 %>%
  group_by(cnt,rep,intAnno2) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom2.intAnno2 <- stat2_intAnno2.s %>%
  ggplot(aes(x = intAnno2, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title="Cell Composition of Nb5d.CTL_CKO, intAnno2", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x=intAnno2, y = 100*prop, color= cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=2.15,
             stroke=0.15, show.legend=F)

pcom2.intAnno2

sig.test

glm.nb - anova.LRT

Sort.N <- c("Nb5d.CTL","Nb5d.CKO")

glm.nb_anova.2.intAnno1 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat2_intAnno1.s_N <- stat2_intAnno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat2_intAnno1.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat2_intAnno1.s_N$total <- total.N[paste0(stat2_intAnno1.s_N$cnt,
                                            "_",
                                            stat2_intAnno1.s_N$rep),"total"]
      
      glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat2_intAnno1.s_N$intAnno1)){
        glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.2.intAnno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.intAnno1)))
rownames(glm.nb_anova.2.intAnno1_df) <- names(glm.nb_anova.2.intAnno1)
colnames(glm.nb_anova.2.intAnno1_df) <- gsub("X","C",colnames(glm.nb_anova.2.intAnno1_df))
glm.nb_anova.2.intAnno1_df
##                           EMN1      EMN2      EMN3       EMN4      EMN5
## Nb5d.CTLvsNb5d.CKO 0.006542252 0.5432009 0.9819977 0.02373406 0.3750407
##                         IMN1      IMN2      IMN3      IMN4       IN1       IN2
## Nb5d.CTLvsNb5d.CKO 0.4678658 0.9067255 0.6576028 0.7693935 0.7436069 0.3062167
##                          IN3     IPAN1     IPAN2     IPAN3     IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7938345 0.7798518 0.3699458 0.4239918 0.7177276
#options (scipen = 999)
round(glm.nb_anova.2.intAnno1_df,4)
##                      EMN1   EMN2  EMN3   EMN4  EMN5   IMN1   IMN2   IMN3   IMN4
## Nb5d.CTLvsNb5d.CKO 0.0065 0.5432 0.982 0.0237 0.375 0.4679 0.9067 0.6576 0.7694
##                       IN1    IN2    IN3  IPAN1  IPAN2 IPAN3  IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7436 0.3062 0.7938 0.7799 0.3699 0.424 0.7177
Sort.N <- c("Nb5d.CTL","Nb5d.CKO")

glm.nb_anova.2.intAnno2 <- list()

for(x1 in 1:2){
  for(x2 in 1:2){
    if(x2>x1){
      
      stat2_intAnno2.s_N <- stat2_intAnno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
      
      total.N <- stat2_intAnno2.s_N %>%
        group_by(cnt, rep) %>%
        summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
      
      rownames(total.N) <- paste0(as.character(total.N$cnt),
                                  "_",
                                  as.character(total.N$rep))
      
      stat2_intAnno2.s_N$total <- total.N[paste0(stat2_intAnno2.s_N$cnt,
                                            "_",
                                            stat2_intAnno2.s_N$rep),"total"]
      
      glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
      
      for(AA in levels(stat2_intAnno2.s_N$intAnno2)){
        glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
          aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
          aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
          aaa.anova <- anova(aaa.glmnb, test = "LRT")
          aaa.anova$'Pr(>Chi)'[2]
        }, error=function(e){
        
          tryCatch({
            aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
            aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
            aaa.anova <- anova(aaa.glmnb, test = "LRT")
            aaa.anova$'Pr(>Chi)'[2]
          }, error=function(e){
            
            tryCatch({
              aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
              aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
              aaa.anova <- anova(aaa.glmnb, test = "LRT")
              aaa.anova$'Pr(>Chi)'[2]
            }, error = function(e){
              NA
            })     
          })
        })
      }
      glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]]) 
    }
  }
}
glm.nb_anova.2.intAnno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.intAnno2)))
rownames(glm.nb_anova.2.intAnno2_df) <- names(glm.nb_anova.2.intAnno2)
colnames(glm.nb_anova.2.intAnno2_df) <- gsub("X","C",colnames(glm.nb_anova.2.intAnno2_df))
glm.nb_anova.2.intAnno2_df
##                           EMN1      EMN2      EMN3       EMN4      EMN5
## Nb5d.CTLvsNb5d.CKO 0.006542252 0.5432009 0.9819977 0.02373406 0.3750407
##                         IMN1      IMN2      IMN3      IMN4       IN1       IN2
## Nb5d.CTLvsNb5d.CKO 0.4678658 0.9067255 0.6576028 0.7693935 0.7436069 0.3062167
##                          IN3   IPAN1.1   IPAN1.2   IPAN2.1  IPAN2.2     IPAN3
## Nb5d.CTLvsNb5d.CKO 0.7938345 0.6342571 0.9505649 0.1930516 0.668667 0.4239918
##                        IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7177276
#options (scipen = 999)
round(glm.nb_anova.2.intAnno2_df,4)
##                      EMN1   EMN2  EMN3   EMN4  EMN5   IMN1   IMN2   IMN3   IMN4
## Nb5d.CTLvsNb5d.CKO 0.0065 0.5432 0.982 0.0237 0.375 0.4679 0.9067 0.6576 0.7694
##                       IN1    IN2    IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2 IPAN3
## Nb5d.CTLvsNb5d.CKO 0.7436 0.3062 0.7938  0.6343  0.9506  0.1931  0.6687 0.424
##                     IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7177

all markers

intAnno1

# find markers for every cluster compared to all remaining cells
Idents(GEX.seur) <- "intAnno1"

GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")
#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.01,
#                                  assay = "SCT",
#                                  test.use = "MAST",
#                                  logfc.threshold = 0.1)
GEX.markers.pre <- read.table("Baf53cre_Nb.markers.SCT_intAnno1.202402.csv", header = TRUE, sep = ",")
#GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
                          levels = levels(GEX.seur$intAnno1))


markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
ttt = 471
ttt/60
## [1] 7.85
ttt/64
## [1] 7.359375
ttt/65
## [1] 7.246154
pp.t60 <- list()

for(i in 1:8){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t60[(60*i-59):(60*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t60
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

ttt = 840
ttt/60
## [1] 14
ttt/64
## [1] 13.125
ttt/65
## [1] 12.92308
pp.t120 <- list()

for(i in 1:14){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t120[(60*i-59):(60*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t120
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

intAnno2

# find markers for every cluster compared to all remaining cells
Idents(GEX.seur) <- "intAnno2"

#GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")
#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.05,
#                                  assay = "SCT",
#                                  test.use = "MAST",
#                                  logfc.threshold = 0.25)
GEX.markers.pre <- read.table("Baf53cre_Nb.markers.SCT_intAnno2.202402.csv", header = TRUE, sep = ",")
#GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
                          levels = levels(GEX.seur$intAnno2))


markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
ttt = 485
ttt/60
## [1] 8.083333
ttt/64
## [1] 7.578125
ttt/65
## [1] 7.461538
pp.t60 <- list()

for(i in 1:8){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t60[(64*i-63):(64*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t60
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

ttt = 875
ttt/60
## [1] 14.58333
ttt/64
## [1] 13.67188
ttt/65
## [1] 13.46154
pp.t120 <- list()

for(i in 1:14){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t120[(64*i-63):(64*i)]))  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t120
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

DEGs

head(GEX.seur@meta.data)
##                        orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1_1 Nb5d.PBS_INF       3257         1801 0.36843721  0.3991403
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS_INF       1511          997 0.66181337  0.4632694
## AAACCCAGTTTGTTCT-1_1 Nb5d.PBS_INF       2855         1577 0.98073555  0.3152364
## AAACCCATCCTAGCCT-1_1 Nb5d.PBS_INF       2433         1451 0.08220304  0.3699137
## AAACCCATCGAAACAA-1_1 Nb5d.PBS_INF       3129         1656 0.12783637  0.4474273
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS_INF       2201         1294 0.22716947  0.2271695
##                           S.Score     G2M.Score Phase      cnt  rep newAnno
## AAACCCACAAGACGAC-1_1  0.011590405 -0.0004169865     S Nb5d.INF rep4    EMN3
## AAACCCAGTGGGCTCT-1_1 -0.024203070  0.0012414826   G2M Nb5d.PBS rep4   IPAN1
## AAACCCAGTTTGTTCT-1_1 -0.013980476  0.0039329410   G2M Nb5d.INF rep1    EMN3
## AAACCCATCCTAGCCT-1_1 -0.028925620 -0.0132582758    G1 Nb5d.INF rep2    EMN1
## AAACCCATCGAAACAA-1_1 -0.008077289 -0.0028336129    G1 Nb5d.PBS rep3   IPAN4
## AAACCCATCGGTCAGC-1_1 -0.023612751  0.0327239644   G2M Nb5d.PBS rep4    EMN1
##                         sample tissue nCount_SCT nFeature_SCT condition
## AAACCCACAAGACGAC-1_1 Nb5d.INF4  Ileum       2592         1794   INF_CTL
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS4  Ileum       1694          996   PBS_CTL
## AAACCCAGTTTGTTCT-1_1 Nb5d.INF1  Ileum       2495         1576   INF_CTL
## AAACCCATCCTAGCCT-1_1 Nb5d.INF2  Ileum       2324         1451   INF_CTL
## AAACCCATCGAAACAA-1_1 Nb5d.PBS3  Ileum       2552         1646   PBS_CTL
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS4  Ileum       2171         1293   PBS_CTL
##                      seurat_clusters sort_clusters intAnno1 intAnno2
## AAACCCACAAGACGAC-1_1              11            11     EMN2     EMN2
## AAACCCAGTGGGCTCT-1_1              22            22    IPAN1  IPAN1.1
## AAACCCAGTTTGTTCT-1_1              11            11     EMN2     EMN2
## AAACCCATCCTAGCCT-1_1               4             4     EMN1     EMN1
## AAACCCATCGAAACAA-1_1              19            19    IPAN4    IPAN4
## AAACCCATCGGTCAGC-1_1               8             8     EMN1     EMN1
Idents(GEX.seur) <- "cnt"

#GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")
neur.clusters <- grep("EMN|IMN|IN|IPAN",levels(GEX.seur$intAnno2), value = T)
neur.clusters
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"
#
all_new <- neur.clusters
all_new
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"
#
list_new <- list()
list_new[["All"]] <- all_new

list_new[["EMNs"]] <- grep("EMN",all_new,value = T)
list_new[["IMNs"]] <- grep("IMN",all_new,value = T)

list_new[["IPAN1"]] <- grep("IPAN1",all_new,value = T)
list_new[["IPAN2"]] <- grep("IPAN2",all_new,value = T)

#list_new[["INs"]] <- grep("IN",all_new,value = T)
#list_new[["IPANs"]] <- grep("IPAN",all_new,value = T)

for(NN in all_new){
  list_new[[NN]] <- NN
}

names_new <- names(list_new)
list_new
## $All
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"  
## 
## $EMNs
## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5"
## 
## $IMNs
## [1] "IMN1" "IMN2" "IMN3" "IMN4"
## 
## $IPAN1
## [1] "IPAN1.1" "IPAN1.2"
## 
## $IPAN2
## [1] "IPAN2.1" "IPAN2.2"
## 
## $EMN1
## [1] "EMN1"
## 
## $EMN2
## [1] "EMN2"
## 
## $EMN3
## [1] "EMN3"
## 
## $EMN4
## [1] "EMN4"
## 
## $EMN5
## [1] "EMN5"
## 
## $IMN1
## [1] "IMN1"
## 
## $IMN2
## [1] "IMN2"
## 
## $IMN3
## [1] "IMN3"
## 
## $IMN4
## [1] "IMN4"
## 
## $IN1
## [1] "IN1"
## 
## $IN2
## [1] "IN2"
## 
## $IN3
## [1] "IN3"
## 
## $IPAN1.1
## [1] "IPAN1.1"
## 
## $IPAN1.2
## [1] "IPAN1.2"
## 
## $IPAN2.1
## [1] "IPAN2.1"
## 
## $IPAN2.2
## [1] "IPAN2.2"
## 
## $IPAN3
## [1] "IPAN3"
## 
## $IPAN4
## [1] "IPAN4"
names_new
##  [1] "All"     "EMNs"    "IMNs"    "IPAN1"   "IPAN2"   "EMN1"    "EMN2"   
##  [8] "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"    "IMN3"    "IMN4"   
## [15] "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2" "IPAN2.1" "IPAN2.2"
## [22] "IPAN3"   "IPAN4"

PBS_INF

test1.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.PBS","Nb5d.INF"))
test1.seur
## An object of class Seurat 
## 47356 features across 8232 samples within 3 assays 
## Active assay: SCT (20434 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap

run

#test1.DEGs_new
# save DEGs' table
df_test1.DEGs_new <- data.frame()
for(nn in names_new){
  df_test1.DEGs_new <- rbind(df_test1.DEGs_new, data.frame(test1.DEGs_new[[nn]],intAnno2=nn))
}

#write.csv(df_test1.DEGs_new, "Baf53cre_Nb.DEGs.PBSvsINF.intAnno2.csv")

stat

df_test1.DEGs_new$intAnno2 <- factor(as.character(df_test1.DEGs_new$intAnno2),
                                     levels = names_new)
cut.padj = 0.05
cut.log2FC = 0.1   
  
cut.pct1 = 0.05

stat_test1a.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1a.DEGs_new
##    intAnno2  cluster up.DEGs
## 1       All Nb5d.PBS     104
## 2       All Nb5d.INF     128
## 3      EMNs Nb5d.PBS      47
## 4      EMNs Nb5d.INF      56
## 5      IMNs Nb5d.PBS      11
## 6      IMNs Nb5d.INF      14
## 7     IPAN1 Nb5d.PBS      43
## 8     IPAN1 Nb5d.INF     127
## 9     IPAN2 Nb5d.PBS      15
## 10    IPAN2 Nb5d.INF      63
## 11     EMN1 Nb5d.PBS      21
## 12     EMN1 Nb5d.INF      31
## 13     EMN2 Nb5d.PBS       6
## 14     EMN2 Nb5d.INF      14
## 15     EMN3 Nb5d.INF       4
## 16     EMN4 Nb5d.INF       2
## 17     IMN1 Nb5d.PBS       3
## 18     IMN1 Nb5d.INF       5
## 19     IMN2 Nb5d.INF       1
## 20     IMN3 Nb5d.PBS       1
## 21     IMN3 Nb5d.INF       1
## 22     IMN4 Nb5d.PBS       1
## 23     IMN4 Nb5d.INF       1
## 24      IN1 Nb5d.INF       1
## 25  IPAN1.1 Nb5d.PBS      17
## 26  IPAN1.1 Nb5d.INF      71
## 27  IPAN1.2 Nb5d.PBS      22
## 28  IPAN1.2 Nb5d.INF      74
## 29  IPAN2.1 Nb5d.PBS       9
## 30  IPAN2.1 Nb5d.INF      60
## 31  IPAN2.2 Nb5d.INF       7
## 32    IPAN4 Nb5d.INF       5
df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

cut.padj = 0.05
cut.log2FC = 0.2   
  
cut.pct1 = 0.05

stat_test1b.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1b.DEGs_new
##    intAnno2  cluster up.DEGs
## 1       All Nb5d.PBS      73
## 2       All Nb5d.INF      98
## 3      EMNs Nb5d.PBS      41
## 4      EMNs Nb5d.INF      53
## 5      IMNs Nb5d.PBS      10
## 6      IMNs Nb5d.INF      13
## 7     IPAN1 Nb5d.PBS      43
## 8     IPAN1 Nb5d.INF     127
## 9     IPAN2 Nb5d.PBS      15
## 10    IPAN2 Nb5d.INF      63
## 11     EMN1 Nb5d.PBS      18
## 12     EMN1 Nb5d.INF      29
## 13     EMN2 Nb5d.PBS       6
## 14     EMN2 Nb5d.INF      14
## 15     EMN3 Nb5d.INF       4
## 16     EMN4 Nb5d.INF       2
## 17     IMN1 Nb5d.PBS       3
## 18     IMN1 Nb5d.INF       5
## 19     IMN2 Nb5d.INF       1
## 20     IMN3 Nb5d.PBS       1
## 21     IMN3 Nb5d.INF       1
## 22     IMN4 Nb5d.PBS       1
## 23     IMN4 Nb5d.INF       1
## 24      IN1 Nb5d.INF       1
## 25  IPAN1.1 Nb5d.PBS      17
## 26  IPAN1.1 Nb5d.INF      71
## 27  IPAN1.2 Nb5d.PBS      22
## 28  IPAN1.2 Nb5d.INF      74
## 29  IPAN2.1 Nb5d.PBS       9
## 30  IPAN2.1 Nb5d.INF      60
## 31  IPAN2.2 Nb5d.INF       7
## 32    IPAN4 Nb5d.INF       5
df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test1, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

plot

#save(test1.DEGs_new, file = "test1.DEGs_new.Rdata")
pp_test1.DEGs.new <- lapply(list_new, function(x){NA})
#
test1.DEGs_new.combine <- test1.DEGs_new
test1.DEGs_new.combine <- lapply(test1.DEGs_new.combine, function(x){
  x[x$cluster=="Nb5d.PBS","avg_log2FC"] <- -x[x$cluster=="Nb5d.PBS","avg_log2FC"]
  x
})

All

pp_test1.DEGs.new$All <- test1.DEGs_new.combine$All %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="All Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$All

EMNs

pp_test1.DEGs.new$EMNs <- test1.DEGs_new.combine$EMNs %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMNs Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMNs

EMN1

pp_test1.DEGs.new$EMN1 <- test1.DEGs_new.combine$EMN1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMN1 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN1

EMN2

pp_test1.DEGs.new$EMN2 <- test1.DEGs_new.combine$EMN2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMN2 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN2

EMN3

pp_test1.DEGs.new$EMN3 <- test1.DEGs_new.combine$EMN3 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMN3 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN3

EMN4

pp_test1.DEGs.new$EMN4 <- test1.DEGs_new.combine$EMN4 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMN4 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN4

IMNs

pp_test1.DEGs.new$IMNs <- test1.DEGs_new.combine$IMNs %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IMNs Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMNs

IMN1

pp_test1.DEGs.new$IMN1 <- test1.DEGs_new.combine$IMN1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IMN1 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN1

IMN2

pp_test1.DEGs.new$IMN2 <- test1.DEGs_new.combine$IMN2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IMN2 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN2

IMN3

pp_test1.DEGs.new$IMN3 <- test1.DEGs_new.combine$IMN3 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IMN3 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN3

IMN4

pp_test1.DEGs.new$IMN4 <- test1.DEGs_new.combine$IMN4 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IMN4 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN4

IN1

pp_test1.DEGs.new$IN1 <- test1.DEGs_new.combine$IN1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IN1 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IN1

IPAN1

pp_test1.DEGs.new$IPAN1 <- test1.DEGs_new.combine$IPAN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1

IPAN1.1

pp_test1.DEGs.new$IPAN1.1 <- test1.DEGs_new.combine$IPAN1.1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.1 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1.1

IPAN1.2

pp_test1.DEGs.new$IPAN1.2 <- test1.DEGs_new.combine$IPAN1.2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.2 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1.2

IPAN2

pp_test1.DEGs.new$IPAN2 <- test1.DEGs_new.combine$IPAN2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2

IPAN2.1

pp_test1.DEGs.new$IPAN2.1 <- test1.DEGs_new.combine$IPAN2.1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.1 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2.1

IPAN2.2

pp_test1.DEGs.new$IPAN2.2 <- test1.DEGs_new.combine$IPAN2.2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.2 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2.2

IPAN4

pp_test1.DEGs.new$IPAN4 <- test1.DEGs_new.combine$IPAN4 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
                               "Nb5d.INF"=color.test1[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN4 Nb5d.PBS vs Nb5d.INF") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN4

neur.ol <- c("All","EMNs","EMN1","EMN2","IMNs","IMN1","IPAN1","IPAN1.1","IPAN1.2","IPAN2","IPAN2.1")
neur.ol
##  [1] "All"     "EMNs"    "EMN1"    "EMN2"    "IMNs"    "IMN1"    "IPAN1"  
##  [8] "IPAN1.1" "IPAN1.2" "IPAN2"   "IPAN2.1"

CTL_CKO

test2.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.CTL","Nb5d.CKO"))
test2.seur
## An object of class Seurat 
## 47356 features across 11186 samples within 3 assays 
## Active assay: SCT (20434 features, 0 variable features)
##  2 other assays present: RNA, integrated
##  3 dimensional reductions calculated: pca, tsne, umap

run

#test2.DEGs_new
# save DEGs' table
df_test2.DEGs_new <- data.frame()
for(nn in names_new){
  df_test2.DEGs_new <- rbind(df_test2.DEGs_new, data.frame(test2.DEGs_new[[nn]],intAnno2=nn))
}

#write.csv(df_test2.DEGs_new, "Baf53cre_Nb.DEGs.CTLvsCKO.intAnno2.csv")

stat

df_test2.DEGs_new$intAnno2 <- factor(as.character(df_test2.DEGs_new$intAnno2),
                                     levels = names_new)
cut.padj = 0.05
cut.log2FC = 0.1   
  
cut.pct1 = 0.05

stat_test2a.DEGs_new <- df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test2a.DEGs_new
##    intAnno2  cluster up.DEGs
## 1       All Nb5d.CTL     159
## 2       All Nb5d.CKO      52
## 3      EMNs Nb5d.CTL      66
## 4      EMNs Nb5d.CKO      20
## 5      IMNs Nb5d.CTL      13
## 6      IMNs Nb5d.CKO       4
## 7     IPAN1 Nb5d.CTL     107
## 8     IPAN1 Nb5d.CKO      30
## 9     IPAN2 Nb5d.CTL      39
## 10    IPAN2 Nb5d.CKO       7
## 11     EMN1 Nb5d.CTL      22
## 12     EMN1 Nb5d.CKO       3
## 13     EMN2 Nb5d.CTL       4
## 14     EMN2 Nb5d.CKO       2
## 15     EMN3 Nb5d.CTL       5
## 16     IMN1 Nb5d.CTL       9
## 17     IMN1 Nb5d.CKO       1
## 18      IN1 Nb5d.CTL       1
## 19  IPAN1.1 Nb5d.CTL      46
## 20  IPAN1.1 Nb5d.CKO      17
## 21  IPAN1.2 Nb5d.CTL      60
## 22  IPAN1.2 Nb5d.CKO      16
## 23  IPAN2.1 Nb5d.CTL      38
## 24  IPAN2.1 Nb5d.CKO       9
## 25  IPAN2.2 Nb5d.CTL       1
df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

cut.padj = 0.05
cut.log2FC = 0.2   
  
cut.pct1 = 0.05

stat_test2b.DEGs_new <- df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat_test2b.DEGs_new
##    intAnno2  cluster up.DEGs
## 1       All Nb5d.CTL      74
## 2       All Nb5d.CKO       8
## 3      EMNs Nb5d.CTL      59
## 4      EMNs Nb5d.CKO       9
## 5      IMNs Nb5d.CTL      13
## 6      IMNs Nb5d.CKO       1
## 7     IPAN1 Nb5d.CTL     107
## 8     IPAN1 Nb5d.CKO      30
## 9     IPAN2 Nb5d.CTL      39
## 10    IPAN2 Nb5d.CKO       7
## 11     EMN1 Nb5d.CTL      21
## 12     EMN1 Nb5d.CKO       2
## 13     EMN2 Nb5d.CTL       4
## 14     EMN2 Nb5d.CKO       1
## 15     EMN3 Nb5d.CTL       5
## 16     IMN1 Nb5d.CTL       9
## 17     IMN1 Nb5d.CKO       1
## 18      IN1 Nb5d.CTL       1
## 19  IPAN1.1 Nb5d.CTL      46
## 20  IPAN1.1 Nb5d.CKO      17
## 21  IPAN1.2 Nb5d.CTL      60
## 22  IPAN1.2 Nb5d.CKO      16
## 23  IPAN2.1 Nb5d.CTL      38
## 24  IPAN2.1 Nb5d.CKO       9
## 25  IPAN2.2 Nb5d.CTL       1
df_test2.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.test2, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))

plot

#save(test2.DEGs_new, file = "test2.DEGs_new.Rdata")
pp_test2.DEGs.new <- lapply(list_new, function(x){NA})
#
test2.DEGs_new.combine <- test2.DEGs_new
test2.DEGs_new.combine <- lapply(test2.DEGs_new.combine, function(x){
  x[x$cluster=="Nb5d.CTL","avg_log2FC"] <- -x[x$cluster=="Nb5d.CTL","avg_log2FC"]
  x
})

All

pp_test2.DEGs.new$All <- test2.DEGs_new.combine$All %>%
  mutate(label=ifelse(((p_val_adj < 1e-5 & avg_log2FC<0) | (p_val_adj < 1e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="All Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$All

EMNs

pp_test2.DEGs.new$EMNs <- test2.DEGs_new.combine$EMNs %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMNs Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMNs

EMN1

pp_test2.DEGs.new$EMN1 <- test2.DEGs_new.combine$EMN1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMN1 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN1

EMN2

pp_test2.DEGs.new$EMN2 <- test2.DEGs_new.combine$EMN2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMN2 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN2

EMN3

pp_test2.DEGs.new$EMN3 <- test2.DEGs_new.combine$EMN3 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="EMN3 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN3

IMNs

pp_test2.DEGs.new$IMNs <- test2.DEGs_new.combine$IMNs %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IMNs Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMNs

IMN1

pp_test2.DEGs.new$IMN1 <- test2.DEGs_new.combine$IMN1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IMN1 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMN1

IN1

pp_test2.DEGs.new$IN1 <- test2.DEGs_new.combine$IN1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IN1 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IN1

IPAN1

pp_test2.DEGs.new$IPAN1 <- test2.DEGs_new.combine$IPAN1 %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1

IPAN1.1

pp_test2.DEGs.new$IPAN1.1 <- test2.DEGs_new.combine$IPAN1.1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.1 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1.1

IPAN1.2

pp_test2.DEGs.new$IPAN1.2 <- test2.DEGs_new.combine$IPAN1.2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN1.2 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1.2

IPAN2

pp_test2.DEGs.new$IPAN2 <- test2.DEGs_new.combine$IPAN2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2

IPAN2.1

pp_test2.DEGs.new$IPAN2.1 <- test2.DEGs_new.combine$IPAN2.1 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.1 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2.1

IPAN2.2

pp_test2.DEGs.new$IPAN2.2 <- test2.DEGs_new.combine$IPAN2.2 %>%
  mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        #!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
                        !grepl("^Rps|^Rpl",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
                               "Nb5d.CKO"=color.test2[2],
                               "None"="grey")) +
  theme_classic() + labs(title="IPAN2.2 Nb5d.CTL vs Nb5d.CKO") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2.2

signature scores

top list

use top120-built top markers

load("I:/Shared_win/projects/20231220_10x_SZJ/analysis_plus/top.markerlist.RData")
data.frame(row.names = "topmarkers",lapply(top.list,length))
##            EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3 IMN4 IN1 IN2 IN3 IPAN1.1
## topmarkers   96   32   37   34   70   71   33   49   53  60  56  63      70
##            IPAN1.2 IPAN2.1 IPAN2.2 IPAN3 IPAN4
## topmarkers      51      43      52    52    50
neur.clusters
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"
for(nn in neur.clusters){
  
  GEX.seur <- add_geneset_score(obj = GEX.seur,
                                Assay = "SCT",
                                geneset = top.list[[nn]],
                                setname = paste0("score.",nn))
  
}
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
FeaturePlot(GEX.seur, features = paste0("score.",neur.clusters), ncol = 5) &
  scale_color_gradientn(colors = rev(mapal))

mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
FeaturePlot(GEX.seur, features = paste0("score.",neur.clusters), ncol = 5, raster = T, pt.size = 3.5) &
  scale_color_gradientn(colors = rev(mapal))

vln.list <- list()

vln.list <- lapply(neur.clusters, function(x){
  
  
  VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.4, 0.4, 0.45),
                               size=2.75
                               )
  
  
})

names(vln.list) <- neur.clusters
cowplot::plot_grid(
  plotlist = vln.list ,
  ncol = 9)

overlapped DEGs

process

cut.padj = 0.05
cut.log2FC = 0.1   
  
cut.pct1 = 0.05

stat_test1a.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(intAnno2,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
#stat_test1a.DEGs_new
head(df_test1.DEGs_new)
##                      p_val avg_log2FC pct.1 pct.2    p_val_adj  cluster
## Ctnna3        8.684419e-51  0.4588045 0.846 0.761 2.129854e-46 Nb5d.PBS
## Malat1        1.335552e-36  0.1410444 1.000 1.000 3.275441e-32 Nb5d.PBS
## AY036118      7.553267e-24  0.2576808 0.788 0.750 1.852439e-19 Nb5d.PBS
## Fgfr2         2.040025e-21  0.2291577 0.838 0.769 5.003161e-17 Nb5d.PBS
## 4930447N08Rik 2.483724e-18  0.2687259 0.726 0.659 6.091333e-14 Nb5d.PBS
## Prkn          2.213112e-17  0.2764906 0.516 0.421 5.427657e-13 Nb5d.PBS
##                        gene intAnno2
## Ctnna3               Ctnna3      All
## Malat1               Malat1      All
## AY036118           AY036118      All
## Fgfr2                 Fgfr2      All
## 4930447N08Rik 4930447N08Rik      All
## Prkn                   Prkn      All
#              
cut.padj = 0.05
cut.log2FC = 0.1   
                
cut.pct1 = 0.05


#
test1_up.list <- lapply(neur.ol,function(x){NA})
names(test1_up.list) <- neur.ol


for(NN in neur.ol){
  
  test1_up.list[[NN]] <- list()
  test1_up.list[[NN]][["Nb5d.PBS"]] <- (df_test1.DEGs_new %>% filter(intAnno2 == NN &
                                                                    cluster == "Nb5d.PBS" &
                                                                    p_val_adj < cut.padj & 
                                                                    abs(avg_log2FC) > cut.log2FC & 
                                                                    pct.1 > cut.pct1))$gene
  test1_up.list[[NN]][["Nb5d.INF"]] <- (df_test1.DEGs_new %>% filter(intAnno2 == NN &
                                                                    cluster == "Nb5d.INF" &
                                                                    p_val_adj < cut.padj & 
                                                                    abs(avg_log2FC) > cut.log2FC & 
                                                                    pct.1 > cut.pct1))$gene
}

#test1_up.list

#
test2_up.list <- lapply(neur.ol,function(x){NA})
names(test2_up.list) <- neur.ol


for(NN in neur.ol){
  
  test2_up.list[[NN]] <- list()
  test2_up.list[[NN]][["Nb5d.CTL"]] <- (df_test2.DEGs_new %>% filter(intAnno2 == NN &
                                                                    cluster == "Nb5d.CTL" &
                                                                    p_val_adj < cut.padj & 
                                                                    abs(avg_log2FC) > cut.log2FC & 
                                                                    pct.1 > cut.pct1))$gene
  test2_up.list[[NN]][["Nb5d.CKO"]] <- (df_test2.DEGs_new %>% filter(intAnno2 == NN &
                                                                    cluster == "Nb5d.CKO" &
                                                                    p_val_adj < cut.padj & 
                                                                    abs(avg_log2FC) > cut.log2FC & 
                                                                    pct.1 > cut.pct1))$gene
}

#test2_up.list
library(UpSetR)
#
xx=neur.ol[1]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[2]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[3]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[4]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[5]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[6]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[7]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[8]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[9]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[10]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

#
xx=neur.ol[11]
  upset(fromList(c(test1_up.list[[xx]],
                   test2_up.list[[xx]])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)  
  grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))

save

test3_up.list <- list()
test3_up.list[["INFxCTL"]] <- list()
test3_up.list[["PBSxCKO"]] <- list()
for(NN in c("All","IPAN1","IPAN1.1","IPAN1.2","IPAN2","IPAN2.1")){
  
  test3_up.list[["INFxCTL"]][[NN]] <- intersect(test1_up.list[[NN]][["Nb5d.INF"]],test2_up.list[[NN]][["Nb5d.CTL"]])
  test3_up.list[["PBSxCKO"]][[NN]] <- intersect(test1_up.list[[NN]][["Nb5d.PBS"]],test2_up.list[[NN]][["Nb5d.CKO"]])
}
test3_up.list
## $INFxCTL
## $INFxCTL$All
##  [1] "Fam155a"       "4930432L08Rik" "Nmu"           "Parm1"        
##  [5] "Cacnb4"        "Pcsk2"         "Calcb"         "Rab27b"       
##  [9] "Nell1"         "Gm30613"       "Tmcc3"         "Gm21847"      
## [13] "Gng2"          "Hnrnpll"       "Cyyr1"         "Grp"          
## 
## $INFxCTL$IPAN1
##  [1] "Parm1"         "Nmu"           "Cadm2"         "Sgip1"        
##  [5] "App"           "Fam155a"       "Dgki"          "4930432L08Rik"
##  [9] "Kcnh8"         "Gng2"          "Cacnb4"        "Hnrnpll"      
## [13] "Col24a1"       "Gm26871"       "Alk"           "Calcb"        
## [17] "Asxl3"         "Nell1"         "Pcsk2"         "Abi3bp"       
## [21] "Igf1r"         "Lingo2"        "Dpyd"          "Tll1"         
## [25] "Negr1"         "Dysf"          "Syt9"          "Gm30094"      
## [29] "Scn3a"         "Epha7"         "Itih5"         "Syt1"         
## [33] "Dgkg"          "Ppm1h"         "Ctnnd2"        "Galnt13"      
## [37] "Gm12709"       "1700111E14Rik" "Epha6"         "Rcan3"        
## [41] "Ptchd4"        "Rab27b"        "Gm38405"       "Ppp1r12b"     
## [45] "Kif5a"         "Pak7"          "Dner"          "Tmcc3"        
## [49] "Fhl1"          "Dclk1"         "Grp"           "Gpc6"         
## [53] "Galm"          "Map3k3"        "Gm21847"       "Rock1"        
## [57] "Kcna4"         "Gm16158"       "Slco3a1"       "Gm49226"      
## [61] "Rad51b"        "Gucy1a2"       "Slco4c1"       "Klf6"         
## [65] "Nell1os"       "Cpne7"         "Gm43391"       "Pcsk2os2"     
## [69] "Gm10791"       "Adk"           "Gm48283"       "Tmem8b"       
## [73] "Unc13b"        "Gm21798"       "Syt14"         "Mapk8"        
## [77] "Npr2"         
## 
## $INFxCTL$IPAN1.1
##  [1] "Nmu"           "Parm1"         "Cadm2"         "Fam155a"      
##  [5] "Sgip1"         "App"           "Kcnh8"         "Dgki"         
##  [9] "4930432L08Rik" "Gng2"          "Hnrnpll"       "Asxl3"        
## [13] "Dysf"          "Dgkg"          "Abi3bp"        "Pcsk2"        
## [17] "Cacnb4"        "Calcb"         "Epha6"         "Lingo2"       
## [21] "Epha7"         "Gm30094"       "Igf1r"         "Ppp1r12b"     
## [25] "Scn3a"         "Ppm1h"         "Nell1"         "Rcan3"        
## [29] "Fhl1"          "Tmcc3"         "Gm12709"       "1700111E14Rik"
## [33] "Grp"           "Rock1"         "Cpne7"         "Gucy1a2"      
## [37] "Klf6"         
## 
## $INFxCTL$IPAN1.2
##  [1] "Parm1"         "Nmu"           "Col24a1"       "App"          
##  [5] "4930432L08Rik" "Alk"           "Sgip1"         "Cadm2"        
##  [9] "Fam155a"       "Dgki"          "Cacnb4"        "Gm26871"      
## [13] "Gng2"          "Nell1"         "Lingo2"        "Hnrnpll"      
## [17] "Negr1"         "Tll1"          "Calcb"         "Rab27b"       
## [21] "Kcnh8"         "Gm38405"       "Igf1r"         "Syt9"         
## [25] "Galnt13"       "Gpc6"          "Dpyd"          "Asxl3"        
## [29] "Pcsk2"         "Itih5"         "Scn3a"         "Kif5a"        
## [33] "Dclk1"         "Abi3bp"        "Ppm1h"         "Pak7"         
## [37] "Gm30094"       "Gm12709"       "1700111E14Rik" "Ctnnd2"       
## [41] "Rad51b"        "Epha7"         "Ptchd4"        "Kcnt2"        
## [45] "Rcan3"         "Tmcc3"         "Epha6"        
## 
## $INFxCTL$IPAN2
##  [1] "Fam155a" "Large1"  "Gm30613" "Nrxn3"   "Pcsk2"   "Gm21847" "Ptprz1" 
##  [8] "Kcnd2"   "Alk"     "Calcb"   "Cacnb4"  "Ppm1h"   "Tmcc3"   "Gm21798"
## [15] "Trpc3"   "Arid5b"  "Nell1"   "Oxr1"    "Galnt13" "Gm30094" "Adcy8"  
## [22] "Kcnq5"   "Rgs7"    "Necab1"  "Dner"    "Dlc1"    "Hnrnpll" "Dmd"    
## [29] "Antxr2"  "Pkp4"   
## 
## $INFxCTL$IPAN2.1
##  [1] "Fam155a" "Large1"  "Gm30613" "Nrxn3"   "Pcsk2"   "Gm21847" "Cacnb4" 
##  [8] "Calcb"   "Gm21798" "Rab27b"  "Alk"     "Ppm1h"   "Ptprz1"  "Nell1"  
## [15] "Tmcc3"   "Arid5b"  "Oxr1"    "Trpc3"   "Gm30094" "Adcy8"   "Kcnq5"  
## [22] "Dclk1"   "Galnt13" "Rgs7"    "Vwc2l"   "Dlc1"    "Necab1"  "Hnrnpll"
## [29] "Antxr2"  "Dner"    "Setbp1" 
## 
## 
## $PBSxCKO
## $PBSxCKO$All
## [1] "Cdh6"    "Agbl4"   "Pdzrn4"  "Rsrp1"   "Zfp804a"
## 
## $PBSxCKO$IPAN1
##  [1] "Zfp804a" "Fgf13"   "Cdh6"    "Tafa1"   "Luzp2"   "Ppp3ca"  "Efr3a"  
##  [8] "Lsamp"   "Arhgap6" "Ano2"    "Rab3c"   "Rbfox1"  "Pdzrn4"  "Cntn5"  
## [15] "Kcnb2"   "Nme6"   
## 
## $PBSxCKO$IPAN1.1
## [1] "Zfp804a" "Fgf13"   "Tafa1"   "Lsamp"   "Arhgap6"
## 
## $PBSxCKO$IPAN1.2
## [1] "Cntn5"   "Zfp804a" "Kcnb2"   "Fgf13"   "Ppp3ca"  "Luzp2"   "Cdh6"   
## [8] "Tafa1"   "Pdzrn4" 
## 
## $PBSxCKO$IPAN2
## [1] "Dgkb"
## 
## $PBSxCKO$IPAN2.1
## [1] "Dgkb" "Cdh6"

plot

length(test3_up.list$INFxCTL$IPAN1)
## [1] 77
vln.DEGs_INFxCTL <- list()
# violin for all DEGs in INFxCTL_IPAN1
#     IPAN1 only
vln.DEGs_INFxCTL[["IPAN1"]] <- lapply(test3_up.list$INFxCTL$IPAN1,function(xx){
  
  VlnPlot(subset(GEX.seur,subset=intAnno1 %in% c("IPAN1")), features = xx, group.by = "cnt", cols = color.cnt, pt.size = 0) +
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1, 1, 1.5),
                               size=2.75
                               )
  
})

names(vln.DEGs_INFxCTL[["IPAN1"]]) <- test3_up.list$INFxCTL$IPAN1
#vln.DEGs_INFxCTL[["IPAN1"]]
cowplot::plot_grid(
plotlist = vln.DEGs_INFxCTL[["IPAN1"]],
ncol = 4
)

length(test3_up.list$INFxCTL$IPAN2)
## [1] 30
# violin for all DEGs in INFxCTL_IPAN2
#     IPAN2 only
vln.DEGs_INFxCTL[["IPAN2"]] <- lapply(test3_up.list$INFxCTL$IPAN2,function(xx){
  
  VlnPlot(subset(GEX.seur,subset=intAnno1 %in% c("IPAN2")), features = xx, group.by = "cnt", cols = color.cnt, pt.size = 0) +
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1, 1, 1.5),
                               size=2.75
                               )
  
})

names(vln.DEGs_INFxCTL[["IPAN2"]]) <- test3_up.list$INFxCTL$IPAN2


# mod one without pvalue
vln.DEGs_INFxCTL[["IPAN2"]][["Fam155a"]] <-  VlnPlot(subset(GEX.seur,subset=intAnno1 %in% c("IPAN2")), features = "Fam155a", group.by = "cnt", cols = color.cnt, pt.size = 0) +
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(3.2, 3.2, 3.8),
                               size=2.75
                               )
#vln.DEGs_INFxCTL[["IPAN2"]]
cowplot::plot_grid(
plotlist = vln.DEGs_INFxCTL[["IPAN2"]],
ncol = 4
)

sig.score

# scoring INFxCTL in IPAN1 and IPAN2
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test3_up.list$INFxCTL$IPAN1,
                              setname = "score.INFxCTL_IPAN1")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test3_up.list$INFxCTL$IPAN2,
                              setname = "score.INFxCTL_IPAN2")
## Summarizing data
vln.INFxCTL <- list()

vln.INFxCTL <- lapply(c("INFxCTL_IPAN1","INFxCTL_IPAN2"), function(x){
  
  
  VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.4, 0.4, 0.45),
                               size=2.75
                               )
  
  
})

names(vln.INFxCTL) <- c("INFxCTL_IPAN1","INFxCTL_IPAN2")
# All
cowplot::plot_grid(
  plotlist = vln.INFxCTL ,
  ncol = 2)

# IPAN1 only
x = "INFxCTL_IPAN1"
vln.INFxCTL_IPAN1 <- VlnPlot(subset(GEX.seur, subset=intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.15,0.75)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.6, 0.6, 0.65),
                               size=2.75
                               )
vln.INFxCTL_IPAN1

# IPAN1 only
x = "INFxCTL_IPAN2"
vln.INFxCTL_IPAN2 <- VlnPlot(subset(GEX.seur, subset=intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.2,1.05)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.85, 0.85, 0.95),
                               size=2.75
                               )
vln.INFxCTL_IPAN2

# scoring DEGs.All
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test1_up.list$All$Nb5d.PBS,
                              setname = "score.All_PBSup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test1_up.list$All$Nb5d.INF,
                              setname = "score.All_INFup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test2_up.list$All$Nb5d.CTL,
                              setname = "score.All_CTLup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test2_up.list$All$Nb5d.CKO,
                              setname = "score.All_CKOup")
## Summarizing data
# scoring DEGs.IPAN1
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test1_up.list$IPAN1$Nb5d.PBS,
                              setname = "score.IPAN1_PBSup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test1_up.list$IPAN1$Nb5d.INF,
                              setname = "score.IPAN1_INFup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test2_up.list$IPAN1$Nb5d.CTL,
                              setname = "score.IPAN1_CTLup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test2_up.list$IPAN1$Nb5d.CKO,
                              setname = "score.IPAN1_CKOup")
## Summarizing data
# scoring DEGs.IPAN2
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test1_up.list$IPAN2$Nb5d.PBS,
                              setname = "score.IPAN2_PBSup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test1_up.list$IPAN2$Nb5d.INF,
                              setname = "score.IPAN2_INFup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test2_up.list$IPAN2$Nb5d.CTL,
                              setname = "score.IPAN2_CTLup")
## Summarizing data
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = test2_up.list$IPAN2$Nb5d.CKO,
                              setname = "score.IPAN2_CKOup")
## Summarizing data
# All_PBSup
x = "All_PBSup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.2),
                               size=2.75
                               )

# All_INFup
x = "All_INFup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.2),
                               size=2.75
                               )

# All_CTLup
x = "All_CTLup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.2),
                               size=2.75
                               )

# All_CKOup
x = "All_CKOup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.2),
                               size=2.75
                               )

vln.DEGs_IPAN1 <- list()

# IPAN1_PBSup
x = "IPAN1_PBSup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.15)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1, 1, 0.85),
                               size=2.75
                               )

# IPAN1_INFup
x = "IPAN1_INFup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,0.75)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.55, 0.55, 0.6),
                               size=2.75
                               )

# IPAN1_CTLup
x = "IPAN1_CTLup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.75)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.55, 0.55, 0.58),
                               size=2.75
                               )



# IPAN1_CKOup
x = "IPAN1_CKOup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.45)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.35, 1.35, 1.15),
                               size=2.75
                               )
vln.DEGs_IPAN1
## $IPAN1_PBSup

## 
## $IPAN1_INFup

## 
## $IPAN1_CTLup

## 
## $IPAN1_CKOup

vln.DEGs_IPAN1.1 <- list()

# IPAN1_PBSup
x = "IPAN1_PBSup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.15)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1, 1, 0.85),
                               size=2.75
                               )

# IPAN1_INFup
x = "IPAN1_INFup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,0.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.55, 0.55, 0.6),
                               size=2.75
                               )

# IPAN1_CTLup
x = "IPAN1_CTLup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.55, 0.55, 0.58),
                               size=2.75
                               )



# IPAN1_CKOup
x = "IPAN1_CKOup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.5)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.35, 1.35, 1.15),
                               size=2.75
                               )
vln.DEGs_IPAN1.1
## $IPAN1_PBSup

## 
## $IPAN1_INFup

## 
## $IPAN1_CTLup

## 
## $IPAN1_CKOup

vln.DEGs_IPAN1.2 <- list()

# IPAN1_PBSup
x = "IPAN1_PBSup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.15,1.15)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1, 1, 0.85),
                               size=2.75
                               )

# IPAN1_INFup
x = "IPAN1_INFup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,0.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.55, 0.55, 0.6),
                               size=2.75
                               )

# IPAN1_CTLup
x = "IPAN1_CTLup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.55, 0.55, 0.58),
                               size=2.75
                               )



# IPAN1_CKOup
x = "IPAN1_CKOup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.2,1.5)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.35, 1.35, 1.15),
                               size=2.75
                               )
vln.DEGs_IPAN1.2
## $IPAN1_PBSup

## 
## $IPAN1_INFup

## 
## $IPAN1_CTLup

## 
## $IPAN1_CKOup

vln.DEGs_IPAN2 <- list()

# IPAN2_PBSup
x = "IPAN2_PBSup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.25)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.1, 1, 0.85),
                               size=2.75
                               )

# IPAN2_INFup
x = "IPAN2_INFup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.7, 0.65, 0.75),
                               size=2.75
                               )

# IPAN2_CTLup
x = "IPAN2_CTLup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.72, 0.72, 0.78),
                               size=2.75
                               )



# IPAN2_CKOup
x = "IPAN2_CKOup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.85)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.65, 1.65, 1.55),
                               size=2.75
                               )
vln.DEGs_IPAN2
## $IPAN2_PBSup

## 
## $IPAN2_INFup

## 
## $IPAN2_CTLup

## 
## $IPAN2_CKOup

vln.DEGs_IPAN2.1 <- list()

# IPAN2_PBSup
x = "IPAN2_PBSup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.25)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.1, 1, 0.85),
                               size=2.75
                               )

# IPAN2_INFup
x = "IPAN2_INFup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.7, 0.65, 0.75),
                               size=2.75
                               )

# IPAN2_CTLup
x = "IPAN2_CTLup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.72, 0.72, 0.78),
                               size=2.75
                               )



# IPAN2_CKOup
x = "IPAN2_CKOup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.85)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.65, 1.65, 1.55),
                               size=2.75
                               )
vln.DEGs_IPAN2.1
## $IPAN2_PBSup

## 
## $IPAN2_INFup

## 
## $IPAN2_CTLup

## 
## $IPAN2_CKOup

vln.DEGs_IPAN2.2 <- list()
  
# IPAN2_PBSup
x = "IPAN2_PBSup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.1)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.95, 0.85, 0.75),
                               size=2.75
                               )

# IPAN2_INFup
x = "IPAN2_INFup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.5, 0.54, 0.58),
                               size=2.75
                               )

# IPAN2_CTLup
x = "IPAN2_CTLup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.7)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.5, 0.45, 0.58),
                               size=2.75
                               )



# IPAN2_CKOup
x = "IPAN2_CKOup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.65)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.45, 1.45, 1.35),
                               size=2.75
                               )
vln.DEGs_IPAN2.2
## $IPAN2_PBSup

## 
## $IPAN2_INFup

## 
## $IPAN2_CTLup

## 
## $IPAN2_CKOup

check IEGs

IEGs.1 <- c("Fos","Jun","Egr1","Arc",
          "Npas4","Homer1","Nptx2","Crem",
          "Syne1", # CPG2
          "Ptgs2", # COX-2
          "Nrn1", # CPG15
          #"Nptx2", # Narp
          "Bdnf", 
          "Plk2", # SNK
          "Plat"  # tPA
          )

IEGs.1 <- IEGs.1[IEGs.1 %in% rownames(GEX.seur)]
IEGs.1
##  [1] "Fos"    "Jun"    "Egr1"   "Arc"    "Npas4"  "Homer1" "Crem"   "Syne1" 
##  [9] "Nrn1"   "Bdnf"   "Plk2"   "Plat"
length(IEGs.1)
## [1] 12
sum(IEGs.1 %in% rownames(GEX.seur))
## [1] 12
pm.All_new.IEGs <- DotPlot(GEX.seur, features = IEGs.1, group.by = "intAnno2",
        cols = c("midnightblue","darkorange1")) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno2")
pm.All_new.IEGs

# scoring IEGs
GEX.seur <- add_geneset_score(obj = GEX.seur,
                              Assay = "SCT",
                              geneset = IEGs.1,
                              setname = "score.IEGs")
## Summarizing data
neur.clusters
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"
 VlnPlot(GEX.seur, features = paste0("score.","IEGs"), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "All neurons") &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.15,0.3)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.15, 0.15, 0.2),
                               size=2.75
                               )
## Warning: Removed 1 rows containing missing values (geom_point).

list_new
## $All
##  [1] "EMN1"    "EMN2"    "EMN3"    "EMN4"    "EMN5"    "IMN1"    "IMN2"   
##  [8] "IMN3"    "IMN4"    "IN1"     "IN2"     "IN3"     "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3"   "IPAN4"  
## 
## $EMNs
## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5"
## 
## $IMNs
## [1] "IMN1" "IMN2" "IMN3" "IMN4"
## 
## $IPAN1
## [1] "IPAN1.1" "IPAN1.2"
## 
## $IPAN2
## [1] "IPAN2.1" "IPAN2.2"
## 
## $EMN1
## [1] "EMN1"
## 
## $EMN2
## [1] "EMN2"
## 
## $EMN3
## [1] "EMN3"
## 
## $EMN4
## [1] "EMN4"
## 
## $EMN5
## [1] "EMN5"
## 
## $IMN1
## [1] "IMN1"
## 
## $IMN2
## [1] "IMN2"
## 
## $IMN3
## [1] "IMN3"
## 
## $IMN4
## [1] "IMN4"
## 
## $IN1
## [1] "IN1"
## 
## $IN2
## [1] "IN2"
## 
## $IN3
## [1] "IN3"
## 
## $IPAN1.1
## [1] "IPAN1.1"
## 
## $IPAN1.2
## [1] "IPAN1.2"
## 
## $IPAN2.1
## [1] "IPAN2.1"
## 
## $IPAN2.2
## [1] "IPAN2.2"
## 
## $IPAN3
## [1] "IPAN3"
## 
## $IPAN4
## [1] "IPAN4"
length(list_new)
## [1] 23
vln.IEGs <- list()

for(NN in names(list_new)){
vln.IEGs[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = paste0("score.","IEGs"), group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.15,0.3)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.15, 0.15, 0.2),
                               size=2.75
                               )
  
}

#vln.IEGs
cowplot::plot_grid(
  plotlist = vln.IEGs,
  ncol = 6
)

vln.Syne1 <- list()

for(NN in names(list_new)){
vln.Syne1[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Syne1", group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2.5)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(1.25, 1.25, 1.6),
                               size=2.75
                               )
  
}

#vln.Syne1
cowplot::plot_grid(
  plotlist = vln.Syne1,
  ncol = 6
)

vln.Fos <- list()

for(NN in names(list_new)){
vln.Fos[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Fos", group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.5),
                               size=2.75
                               )
  
}

#vln.Fos
vln.Jun <- list()

for(NN in names(list_new)){
vln.Jun[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Jun", group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.5),
                               size=2.75
                               )
  
}

#vln.Jun
vln.Nrn1 <- list()

for(NN in names(list_new)){
vln.Nrn1[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Nrn1", group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.5),
                               size=2.75
                               )
  
}

#vln.Nrn1
vln.Bdnf <- list()

for(NN in names(list_new)){
vln.Bdnf[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Bdnf", group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.25, 0.25, 0.5),
                               size=2.75
                               )
  
}

#vln.Bdnf
vln.Il13ra1 <- list()

for(NN in names(list_new)){
vln.Il13ra1[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Il13ra1", group.by = "cnt", cols = color.cnt, pt.size = 0) +
    
  NoLegend() &
  geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
    stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) &  
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
                                                  c("Nb5d.CTL","Nb5d.CKO"),
                                                  c("Nb5d.INF","Nb5d.CTL")),
                               label.y = c(0.85, 0.85, 0.95),
                               size=2.75
                               )
  
}

#vln.Il13ra1
#saveRDS(GEX.seur, "./sn10x_WYS.sct_anno.s.rds")
#GEX.seur <- readRDS("./sn10x_WYS.sct_anno.s.rds")